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Preface 


This  study  was  conducted  as  part  of  the  Houston-Galveston  Navigation 
Channels,  Texas  Project.  The  three-dimensional  hydrodynamic  and  salinity 
study  presented  herein  was  conducted  by  the  Hydraulics  Laboratory  (HL), 
U.S.  Army  Engineer  Waterways  Experiment  Station  (WES).  The  study  was 
authorized  by  the  U.S.  Army  Engineer  District,  Galveston,  in  conjunction 
with  the  study  cost-sharing  partner,  the  Port  of  Houston,  Texas,  on  3  April 
1990. 

The  study  was  conducted  by  HL  personnel  under  the  general  direction  of 
Messrs.  Frank  A.  Herrmann,  Jr.,  Director,  HL;  Richard  A.  Sager,  Assistant 
Director,  HL;  William  H.  McAnally,  Jr.,  Chief,  Estuaries  Division,  HL;  and 
David  R.  Richards,  Chief,  Estuarine  Simulation  Branch,  Estuaries  Division. 
Mr.  William  D.  Martin,  Chief,  Estuarine  Engineering  Branch,  Estuaries 
Division,  was  the  Program  Manager.  Principal  investigator  and  author  of  this 
report  was  Dr.  Bernard  B.  Hsieh,  Estuarine  Simulation  Branch.  Dr.  Robert 
T.  McAdory,  Estuarine  Engineering  Branch,  provided  three-dimensional 
hydrodynamic  simulation  results. 

At  the  time  of  publication  of  this  report.  Director  of  WES  was  Dr.  Robert 
W.  Whalin.  Commander  was  COL  Bruce  K.  Howard,  EN. 
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1  Introduction 


Background 

Galveston  Bay,  the  largest  estuary  on  the  Texas  coast,  located  in  south¬ 
eastern  Texas  along  the  Gulf  of  Mexico,  is  a  biologically  productive  and 
economically  important  estuary.  The  important  commercial  and  recreational 
fisheries  include  oyster,  shrimp,  crab,  and  various  finfish. 

The  study  of  circulation  and  salinities  in  the  Galveston  Bay  system  is  com¬ 
plex.  A  number  of  physical  processes  operate  in  the  bay,  and  their  relative 
importance  can  vary  both  spatially  and  temporally.  Bathymetry  and  geometry 
of  the  bay,  astronomical  tide-induced  currents,  wind-induced  circulation, 
density  variations  and  resulting  gravitational-induced  currents,  and  freshwater 
inflow  are  major  factors  determining  bay-wide  circulation  and  salinity  pat¬ 
terns.  In  addition,  the  proposed  deepening  and  widening  of  the  Galveston 
Entrance  Channel  and  the  Houston  Ship  Channel  could  affect  both  circulation 
and  salinities  throughout  the  bay  system. 

Houston  Ship  Channel,  the  major  navigation  channel  in  Galveston  Bay, 
transects  the  bay  from  Bolivar  Roads  at  the  entrance  to  Galveston  Bay 
northward  to  Morgans  Point.  The  channel  then  continues  up  to  the  Main 
Turning  Basin  near  the  city  of  Houston.  The  Galveston  Channel  is  much 
shorter  in  length  and  bifurcates  from  the  Houston  Ship  Channel  in  the  Bolivar 
Roads  area.  The  channel  reach  from  the  inlet  to  the  Gulf  of  Mexico  is  known 
as  the  Galveston  Entrance  Channel.  At  present  the  width  of  the  Houston  Ship 
Channel  is  122  m  and  the  depth  of  12  m  (40  ft)  at  mean  low  tide  is  main¬ 
tained  along  most  of  the  route. 

The  U.S.  Army  EngiftfeeL.'W^^^terwiys  Experiment  Station  Hydraulics 
Laboratory  is  responsible  for  numerical  mod-A  testing  and  long-  and  short¬ 
term  field  data  collection  to  assess  the  potential  etfectO' of-the  channel -enlarge¬ 
ment  on  the  hydrodynamics  and  transport  of  the  system.  Berger  et  al.  (in 
preparation),  Fagerburg  et  al.  (1994),  and  Lin  (1992)  discuss  other  aspects  of 

this  project. 

With  recent  advances  in  computer  technology,  long-term  (1-year)  simula¬ 
tion  of  complicated  estuaries  with  large-scale  multidimensional  numerical 
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models  has  become  possible.  However,  these  numerical  solutions  still  require 
the  use  of  substantial  supercomputer  central  processing  unit  time.  An  alter¬ 
native  method,  capable  of  analyzing  changes  at  individual  points  in  an  estuary, 
as  opposed  to  the  global  solutions  generated  by  the  numerical  models,  is 
desirable.  This  alternative  method  must  also  be  able  to  be  run  in  a  short 
period  of  time.  A  system  response  approach  using  the  input/output  relation¬ 
ships  from  a  numerical  model  to  describe  the  dynamic  behavior  of  tidal  hydro- 
dynamic  phenomena  can  be  used  to  play  this  role.  The  system  response  func¬ 
tions  from  a  numerical  model  verified  for  a  particular  location  can  be  used  to 
simulate  the  resulting  output  function,  such  as  change  in  salinity,  when  input 
forcing  functions,  such  as  tidal  variation  and  freshwater  inflow,  change. 

This  approach  was  used  to  address  the  salinity  response  due  to  freshwater 
inflow  changes  for  16  selected  locations  in  Galveston  Bay,  Texas  (Figure  1). 
The  system  model  base  was  constructed  by  selecting  node  points  from  three- 
dimensional  (3-D)  numerical  hydrodynamic  model  results.  The  annual 
numerical  simulation  of  both  base  geometry  (12-m-deep  (40-ft-deep)  channel) 
and  project  conditions  (13.7-m-deep  (45-ft-deep)  channel)  for  1990  medium- 
flow  conditions  was  used  to  construct  the  system  response  function.  Three 
major  tributaries  (Trinity  River,  San  Jacinto  River,  and  Buffalo  Bayou)  were 
considered  as  primary  freshwater  inflow  points. 


Objective 

The  objective  of  this  work  was  to  use  the  capability  of  system  simulation 
techniques  to  evaluate  salinity  changes  at  16  selected  points  in  the  bay  for 
conditions  projected  for  medium  freshwater. inflow  conditions  in  1999. 
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Figure  1 .  Location  of  time  series  points 
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2  System  Design  and 

Modeling  Considerations 


Physical  Considerations 


Motion  in  an  estuary  is  composed  of  physical  processes  that  occur  over  a 
range  of  frequencies.  In  the  past,  researchers  have  decomposed  time- 
dependent  velocities,  tides,  etc.  into  the  frequency  domain  and  derived  linear 
transfer  functions  to  fill  in  missing  data.  However,  interactions  among  a  large 
number  of  forcing  functions  (e.g.,  tides,  winds,  and  freshwater  inflows)  are 
sometimes  highly  nonlinear.  Therefore,  transfer  functions  based  on  nonlinear 
theory  with  multiple  inputs  are  required. 

Physically,  the  offshore  tide  becomes  strongly  distorted  as  it  propagates 
into  a  shallow  area  system.  This  is  the  result  of  finite  amplitude  effects 
caused  by  friction,  nonlinear  advection,  and  interactions  with  channel 
geometry  as  the  tide  oscillates  within  the  estuary .  These  tidal  asymmetries 
cause  falling  and  rising  tides  to  be  unequal  in  duration,  resulting  in  tidal 
curves  and  salinity  variation  within  the  tidal  cycle.  This  distortion  can  be 
represented  as  the  nonlinear  growth  of  the  astronomical  constituents.  In  addi¬ 
tion,  the  variation  in  freshwater  inflow  changes  the  amplitude  and  phase  of 
tidal  fluctuations. 

Since  tidal  energy  distribution  is  concentrated  in  several  very  narrow 
frequency  bands,  the  frequency  domain  approach  can  capture  characteristics 
which  the  time  domain  approach  could  miss.  For  instance,  using  a  1-hr  time 
interval  in  the  tidal-induced  environment,  semidiurnal  components  only  appear 
in  the  wavelength  with  periods  of  12  and  hr  n^ajor  M2  component, 
with  a  period  of  12.42  br  .iui  d,  significantly  accounted  for  in  the 

variation. 


System  Simulation 


The  governing  equations  of  estuary  flow  and  transport  are  usually  solved 
by  numerical  methods.  Before  these  numerical  models  can  be  used  to  address 
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real-world  problems,  they  must  be  verified.  This  requires  that  field  data  taken 
at  interior  points  agree  with  the  computational  results  from  the  model.  The 
system  parameters  then  are  saved  for  addressing  the  management  issues  such 
as  salinity  changes  at  points  within  the  estuary.  This  process  is  well  known  as 
simulation  via  numerical  modeling. 

For  simulation  via  the  system  model,  system  identification  is  based  on  the 
known  input  and  output  series.  During  this  identification  process,  the 
frequency  response  function  (FRF)  links  the  relationships  between  input  series 
and  output  functions.  The  FRF  is  obtained  by  solving  a  system  matrix  after 
converting  all  input/output  series  from  the  time  domain  to  the  frequency 
domain  via  the  FFT.  This  FRF  usually  covers  the  frequency  band  over  the 
entire  simulation  period.  For  example,  a  discrete  one-sided  spectrum  nor¬ 
mally  covers  the  frequencies  from  0.000115  to  0.5  for  a  year-long  record  of 
data  taken  at  hourly  intervals.  The  predicted  output  then  is  produced  in  the 
frequency  domain  and  converted  by  taking  this  FRF  via  the  inverse  FFT  back 
to  the  time  domain.  The  system  model  verification  process  consists  of  adjust¬ 
ing  the  band  width  and  other  system  parameters,  such  as  maximum  lag, 
integration  of  frequency  range,  and  memory  lengths  of  response  function. 
These  parameters  are  related  to  simulation  accuracy  during  the  FFT  integra¬ 
tion  processes.  The  system  simulation  is  performed  by  taking  the  best  fit  FRF 
of  the  system  through  the  proposed  input  series  of  management  concern. 
Processes  for  comparing  these  two  simulations  are  shown  in  Figure  2. 


Moving  Response  with  Fixed  Input  Design 

For  an  estuary  system,  the  salinity  level  (output)  at  one  particular  location 
of  a  computational  domain  is  approximately  a  function  of  surface  elevation  at 
the  tidal  boundary,  freshwater  discharge  at  the  river  boundary,  and  wind  forc¬ 
ing  (inputs)  on  the  water  surface.  Therefore,  this  system  can  be  classified  as 
a  multiple  input  system.  In  order  to  incorporate  the  system  model  with  the 
numerical  model  input/output  structure,  a  moving  response  function  with  fixed 
inputs  designed  for  the  model  was  proposed  as  shown  in  Figure  3.  Thus 
each  system  output  corresponds  to  one  computational  node  from  the  numerical 
model.  The  advantage  of  this  design  is  that  the  location  of  the  output  cell  of 
interest  can  change  but  the  input  functions  remain  the  same.  This  provides  a 
relative  index  to  examine  the  response  "strength"  through  the  computational 
domain  due  to  the  node  chosen. 


Group  Time  Delays  and  Nonlinearity 

In  considering  the  travel  time  between  the  input  and  output  series,  the 
phase  difference  is  called  the  time  delay.  Time  delays  are  dependent  on  forc¬ 
ing  distance  and  physical  structure.  For  solving  a  multiple  input  system,  a 
group  of  time  delays  are  involved.  This  produces  very  large  bias  errors  when 
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observed  Numerical  observed 

inputs  ^  Model  output 
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Figure  3.  Moving  response  function  with  fixed  input  structure 


performing  cross-spectral  estimates  and  inverse  FFT.  A  special  treatment 
other  than  simply  reducing  the  width  of  the  frequency  band  or  conducting 
multiple  time-shifting  processes  is  required.  As  described  before,  nonlinearity 
in  the  system  renders  the  traditional  assumption  of  multiple  linear  systems 
invalid.  Solution  of  these  technical  difficulties  is  briefly  summarized  as 
follows:  (a)  the  difficulties  due  to  group  time  delays  and  nonlinearity  are 
eliminated  by  decomposing  the  input  and  output  series  into  a  series  of  parallel 
subsystem  paths  with  equivalent  computational  schemes,  (b)  a  finite-memory 
nonlinear  system  is  inserted  between  the  input  series  and  subtransfer  function. 
It  implies  that  the  multiple  frequency  response  function  (MFRF)  is  the  sum¬ 
mation  of  each  subfrequency  response  function,  (c)  any  time  delay  of  each 
subsystem  can  be  corrected  individually. 
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System  Model  Develop¬ 
ment,  Mathematical 
Description 


Generally,  in  an  estuary  system,  measurable  real  input  data,  such  as 
boundary  forcings,  are  transformed  into  measurable  or  computed  output  data. 
As  a  rule,  such  a  system  is  a  multiple  one,  with  several  inputs  and  several 
outputs.  In  general,  all  the  input/output  series  are  assumed  as  stationary 
series.  In  the  present  case,  a  system  with  several  inputs  and  one  output  is 
considered  a  so-called  multiple  input/single  output  system.  By  way  of  the 
mathematical  system  model,  measured  input  data  are  transformed  into 
fictitious,  computed  output  data.  Under  the  conditions  of  an  approximately 
time-invariant  system,  the  mathematical  operation  which  can  simulate  the  salt 
transport  within  an  estuary  can  be  represented  in  the  form  of  the  convolution 
integral.  The  output  y(t)  (salinity  level  at  one  particular  node  of  the  numerical 
finite-element  hydrodynamic  model)  is,  therefore. 


yit)  =  f  h^T)  xft-T)  d  T+e{t) 
7=1  J 


The  Xj(t-T)  represents  the  input  series  I  at  time  level  t.  The  hj(T)  are  called 
partial  impulse  response  functions,  which  describe  the  estuary  system 
unambiguously  and  completely.  The  function  e(t)  is  an  uncorrelated  noise 
term  which  arises  because  the  input  and  output  variables  may  not  be  well- 
controlled.  If  this  error  term  is  neglected,  the  impulse  response  functions  can 
be  estimated  from  the  values  of  Xj{t)  and  y(t)  assuming  that  the  data  series  are 
realizations  of  stationary  stochastic  processes. 

For  determining  the  impulse  response  functions  from  the  design  sysiem  in 
the  frequency  domain,  it  is  necessary  to  transform  the  input  and  output  func¬ 
tions  from  the  time  domain  into  the  frequency  domain.  To  simplify  the 
mathematical  representation,  a  system  with  one  input  and  one  output  may  be 
considered: 
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y{t)  =  f  hit)  x{t-T)  dr  +  e{f) 


(2) 


Neglecting  e{t),  and  using 


Xif)  =  |x(0  dt 

Y(f)  =  jyW  dt 


H(f)  =  |/J(0  dt 


as  FFT  ofx(t),  y(t),  and  h(t),  respectively,  with  i  =  \fT ,  the  following 
relation  can  be  developed: 


Y(f)  =  H(f)  ■  X(f) 


Equation  6  replaces  the  convolution  integral  (1)  by  a  simple  algebraic 
equation.  The  FFT  H(f)  of  the  impulse  response  function  h(t)  is  called  the 
frequency  response  function.  Inverse  Fourier  transformation  of  H(f)  again 
results  in  the  impulse  response  function 


hit)  =  I  Hif)  dt 


It  is  noted  that  Equation^  ^  iind  /  are  a  transformation  pair.  It  is 
concluded  that  a  liiicar  iime-invariant  system  consists  of  an  input  series,  output 
series  .ad  response  function  if  either  H(f)  or  h(t)  is  given.  It  can  be  shown 
that  the  frequency  response  function  may  be  estimated  with  the  aid  of  the 
cross  spectrum  (or  cross-spectral  density  function)  between  input  and  output 
and  the  power  spectrum  of  the  input  (auto-spectral  density  function). 


Chapter  3  System  Model  Development,  Mathematical  Description 


where  G^(f)  and  G^(f)  are  the  Fourier  transforms  of  the  estimate  of  the 
crossvariance  function  r^(t)  and  of  the  autocovariance  function 
respectively.  For  discrete  data  series  with  a  finite  length,  the  estimate  of  the 
true  spectrum  at  discrete  frequencies  can  be  computed  according  to  selected 
methods  such  as  covariance  function  and  Fast  Fourier  Transform  (FFT).  The 
associated  ordinary  coherence  function  estimate  is 


r^xy(f)  = 


IGxy{f)f 
Gxx(f).  Gyy(f) 


(9) 


where  Gyy(f)  is  a  smooth  estimate  of  the  output  cross-spectrum  function.  The 
range  of  the  coherence  function  is  between  0  and  1 .  For  an  ideal  transfer 
function,  the  coherence  function  is  unity  over  the  entire  frequency  band. 

The  smoothed  value  of  the  true  cross-spectrum  is  a  complex-valued 
quantity  defined  by  the  co-spectrum  component  and  quad-spectrum 
component  as; 


G:cy(f)  =  C^(/)  -  i  Qxy(f) 


(10) 


The  one-sided  cross  spectrum  may  be  presented  in  complex  polar  notation  as 


G^(f)  =  /  G^(f)I 


(11) 


where  the  absolute  value  and  phase  angle  are  determined  by 
IG^(f)I  =  SQRT  (C^xyif)  +  0  V(/)) 


The  absolute  value  IGxy(f)I  is  called  the  gain  factor  and  the  assQciatcd 
phase  angle  Pxy(f)  is  called  the  phase  factor.  In  these  terms,  the  frequency 
response  function  takes  on  a  direct  physical  interpretation  and  gives  an  insight 
into  the  behavior  of  the  underlying  system.  Especially,  the  phase  factor  gives 
a  useful  measure  for  the  time  delay  through  input-output  series  at  any  given 
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frequency  f.  Similarly,  the  relationships  between  Equations  10  and  13  can  be 
combined  with  Equation  8  (since  G^(x)  is  the  auto  spectrum  and  the  real 
number)  to  estimate  the  system  gain  factor  and  system  phase  factor  of 
frequency  response  function  H(f). 

According  to  Equation  1,  it  is  necessary  to  transform  the  frequency 
response  function  back  into  the  time  domain.  For  a  simple  single-input  and 
single-output  (SISO)  system,  the  discrete  expression  of  Equation  1  can  be 
represented  by  taking  the  inverse  Fourier  transform  from  the  estimate  of 
Equation  9. 


U 

y{t)  =  Y.  m)  ■  xit-k)) 

k=-p 

for 

M+1  <t<n-p 


(14) 


The  output  estimate  y(t)  is  the  system  response  estimation  due  to  input 
series  x(t).  The  procedures  for  computing  Equations  3  through  14  constitute 
the  simplest  linear  system  estimate  of  output  function  y(t)  via  frequency 
domain  approach.  The  important  outputs  from  these  computations  are  the 
frequency  response  function,  the  system  gain  factor,  the  system  phase  factor, 
the  coherence  function,  and  the  output  estimate. 

For  modeling  a  real  physical  condition,  a  more  complicated  system  is 
considered.  Usually,  a  nonlinear  system  with  multiple  inputs  is  used  to 
describe  the  system  behavior.  A  model  of  a  linear  system  responding  to  N 
inputs  is  developed  by  extending  the  SISO  relationship  as 

N  N 

y(f)  =  Y  ^  E 

7=1  7=1 


The  system  of  linear  equations 
Gxy(f)  +  Gxx(J)  .  H(f) 


{Hjif)  Xj  (f)) 


(15) 


can  be  written  in  matrix  notation 


(16) 


where 


H(f)  =  [Hiv. .  HN(f)f 


(17) 
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Gxy(j)  +  [G^^yif), 


(18) 


G  „(/)  G  i2(/) -  G  ^j4f) 

Gxx(f)  -  :  :  '■ 

Gni^  -  ^Aw(/)  (19) 


where 

H(f)  =  N-dimensional  vector  of  the  frequency  response  functions  for  one 
frequency  / 

Gxyif)  =  N-dimensional  vector  of  the  cross  spectrum  of  all  inputs  with  the 
output,  for  one  frequency  /  (output  spectrum  matrix) 

Gxx^  =  NxN  matrix  of  the  cross  and  auto  spectrum  of  all  inputs  with 
each  other,  for  one  frequency  /  (input  spectrum  matrix) 

The  multiple  frequency  response  function  H(f)  is  solved  by 
H{f)  =  G'^(f) .  G^(f) 


For  solving  a  multiple  linear  system,  two  solution  options  exist;  namely, 
the  direct  inverse  matrix  computation  method,  or  the  parallel  SISO  linear 
system  decomposition  method.  While  the  first  method  involves  more  compli¬ 
cated  techniques  for  solving  the  inverse  matrix,  such  as  partitioning  the  com¬ 
plex  variable,  the  second  method  offers  an  easier  approach  to  the  computa¬ 
tions  but  it  requires  removal  of  the  dependency  among  inputs.  Usually, 
additional  transformation  is  required  to  convert  correlated  inputs  to  uncor¬ 
related  inputs.  Under  the  second  approach,  a  new  multiple  linear  system  is 
formed  by  the  multiple  frequency  response  function  L(f),  and  the  matrices  of 
Gm(f)  and  Guy0. 


L(f)  =  [Ll(f),  L2(f), . 

. .  Ln(f)f 

(21) 

Guy(f)  =  [Guiyif), - 

(22) 
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^  ulul^  ^  ulu2^  ^  uluN^ 

Guu(f)  =  ^  /A 

^uNul^  ^uNu2^  ^uNuhfj) 


(23) 


To  solve  this  new  linear  system,  extensive  computation  is  required  to  fill 
both  the  output  spectrum  matrix  and  the  input  spectrum  matrix.  However, 
one  alternative  is  to  solve  this  transform  matrix  Lxx{f)  and  use  computed 
GxcOO  and  Gxy^  to  estimate  the  solution  via  linear  combination 
technique. 


U(f)  =  [Ul(f),  U2(f),....,  Un(f)f 


(24) 


where 


Uj(f)  =  Xj(f)  -  (Lij  Ui  (/))  7  =  1, 

(=1 


N 


(25) 


N 


L(f)  =  Y  (LijHjif))  j  =  1, 


N 


i=l 


(26) 


where  Lii=l  for  any  i  and  Lij  =  0  if7<z 


Lxxfj) 


f-ll(/)  ^12(/) 


(27) 


Lij  =  Gij  .  12...(A^-l)/G,.,..1.2...(iV-l) 
/=1, . N,  7=i,  (i+1), . N 


(28) 
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GjjA...ij-l) 


(29) 


L(/)  =LJf).H{f)=LJf).G-^(f).G^{f) 


(30) 


The  multiple  coherence  function  estimate  for  this  system  is 


r,„2  =  [Li(/).G^„,(/)  + . +  L^.G^^MIGyyif) 


(31) 


where 

Gyi-Gyi-Lj(f).G,ji  /=1, . N 


(32) 


The  multiple  coherence  function  estimates  represent  the  fraction  of  output 
power  at  frequency /which  can  be  attributed  to  linear  relationships  in  the 
system. 

For  a  multiple  nonlinear  system,  a  decomposition  procedure  must  be  con¬ 
ducted.  Bendat  (1990)  has  proven  that  a  combination  of  a  finite-memory 
nonlinear  system  with  a  parallel  linear  system  is  equivalent  to  a  multiple  non¬ 
linear  system.  A  finite-memory  nonlinear  system  could  be  as  simple  as  the 
second  order  of  polynomial  function  between  one  of  the  input  series  and  out¬ 
put  function.  The  mathematical  detail  can  be  found  in  Bendat  (1990). 
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4  System  Model  Construction 


Define  Input/Output  Series 


Using  physical  realization  to  define  the  input/output  series  is  the  first  step 
to  construct  the  system  model.  The  energy  driving  sources  usually  are  used  to 
represent  the  input  series  and  the  resulting  physical  parameters  are  used  to 
prescribe  the  output  series.  The  actual  number  of  input  series  may  need  to  be 
identified  by  computing  the  cross  spectrum  or  correlation  coefficients  for  each 
pair  of  input/output  series.  This  allows  the  model  structure  to  remain  simple 
but  include  all  significant  components.  For  a  strong  cyclic  relationship, 
harmonic  analysis  can  aid  in  understanding  the  significant  components.  It  is 
suggested  that  the  better  representation  is  to  keep  the  model  structure  as 
simple  as  possible. 


Trend  Removal 

Unlike  general  numerical  model  construction,  the  system  model  is  assumed 
to  be  stationary  over  the  entire  time  horizon.  It  implies  that  the  trend  removal 
process  is  applied  to  both  input  and  output  series.  Trends  that  have  been 
removed  are  added  back  during  the  simulation  stage.  The  most  common 
technique  for  trend  removal  is  to  fit  a  low-order  polynomial  to  the  data  using 
the  least  squares  estimation  method.  Although  more  complex  trends  can  be 
removed  by  higher-order  polynomial  fits,  trend  removal  using  an  order  of 
magnitude  greater  than  three  is  generally  not  suggested.  Sometimes,  the 
filtering  of  data  prior  to  more  detailed  analyses  may  be  desired  for  various 
reasons,  including  either  the  isolation  or  elimination  of  periodic  components. 


System  Structure  Design 

The  structure  of  the  system  needs  to  be  identified  by  the  order  of  input 
significance,  the  distribution  of  the  energy  spectrum,  nonlinear  effect,  and  the 
speed  of  reaction  between  input  and  output  series.  Several  candidates  are 
necessary  in  the  initial  stage  of  modeling.  The  nonlinear  transformation  and 
parallel  signal  paths  (Bendat  1990)  need  to  be  defined. 
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Model  Verification 


The  system  model  is  tested  by  determining  the  bandwidth  of  the  frequency 
domain,  maximum  lag  number,  integration  frequency  range,  lag  weight  coef¬ 
ficients,  and  memory  lengths  of  the  response  function.  The  model  is,  there¬ 
fore,  tuned  until  an  acceptable  error  criterion  is  achieved.  The  multiple 
coherence  functions  or  the  time  series  of  output  residuals  are  the  measure  of 
modeling  performance. 


System  Simulation 

Once  the  system  model  has  been  verified,  the  system  coefficients  are  saved 
and  used  to  perform  system  simulation.  These  system  coefficients,  especially 
the  multiple  frequency  response  function,  remain  in  the  system  and  combine 
with  desired  new  input  to  generate  the  simulated  output.  It  needs  to  be  men¬ 
tioned  that  the  driving  mechanisms  of  system  modeling  are  variations  over 
time.  Less  variation  in  the  input  series  will  receive  less  response,  which  will 
show  less  impact  on  the  output  function.  This  consideration  should  be 
included  prior  to  construction  of  the  model. 


Chapter  4  System  Model  Construction 


5  Galveston  Bay  Application 


Galveston  Bay,  as  most  estuary  systems,  receives  boundary  forcings  from 
tidal  propagation,  freshwater  inflows  from  major  tributaries,  and  wind  stress 
on  the  water  surface.  The  3-D  hydrodynamic  numerical  model  was  con¬ 
structed  to  simulate  present  and  future  salinity  conditions  due  to  channel  deep¬ 
ening  alternatives,  wind  stress,  and  freshwater  inflow  conditions.  The  system 
model  used  numerical  results  of  annual  hourly  salinity  from  16  selected  nodes 
from  the  computational  mesh  of  the  3-D  numerical  model  as  output  (Fig¬ 
ure  4).  The  system  model  uses  the  same  types  of  forcing  functions  including 
the  source  tide  on  the  boundary,  x-component  and  y-component  wind  forcing 
and  three  major  tributaries;  namely.  Trinity  River,  San  Jacinto  River,  and 
Buffalo  Bayou  to  build  a  six-input/single-output  (salinity)  nonlinear  system. 

All  numerical  results  for  the  top  layer  and  bottom  layer  at  these  16  locations 
(Figure  4,  except  point  12,  with  only  one  layer)  were  used  to  construct  their 
individual  transfer  function  models  for  both  base  and  project  conditions.  This 
made  62  dynamic  system  models  under  the  1990  and  1999  medium-flow 
conditions.  The  system  simulation  runs  were  performed  after  these  62  models 
were  verified. 

The  first  step  of  the  analysis  was  to  normalize  the  input/output  series  into 
the  same  scale  of  variation.  The  wind  components  and  flow  data  were 
multiplied  by  factors  of  1000.0  and  0.001,  respectively.  Increasing  the  fresh¬ 
water  inflows  causes  decreasing  salinity  concentrations,  which  is  not  consistent 
with  the  fact  that  increasing  the  surface  elevation  at  the  open  (Gulf)  boundary 
causes  increasing  salinity  values.  Therefore,  a  transformation  taking 
reciprocal  values  of  flow  was  performed.  This  was  the  simplest  nonlinear 
transformation  among  a  group  of  candidates.  This  also  assumed  that  the 
system  was  zero-memory  nonlinear  through  each  input/output  series. 

A  long-term  time  series  can  be  decomposed  into  three  major  components  to 
describe  its  dynamic  behavior;  namelv  the  trend  component,  cyclic  or 
periodic  component,  and  stochastiw  component.  Since  the  relationships  among 
the  input  and  output  series  were  constructed  under  a  stationary  process,  the 
trend  had  to  be  remo’  ed  first.  After  several  tests,  a  slight  linear  trend 
component  w.a.:  removed  from  all  series.  The  frequency  domain  approach  can 
simulate*  i.ach  cyclic  component  through  the  prescribed  techniques.  Therefore, 
the  lesiduals  (white  noise)  can  be  estimated  after  determining  both  the  trend 
and  the  periodic  components  through  the  system. 
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Harmonic  analysis  was  used  to  identify  the  rank  of  input  significance  and 
the  distribution  of  the  energy  spectrum.  It  was  found  that  the  diurnal  tidal 
components  contribute  the  most  high-frequency  variation  but  the  subtidal 
component  and  neap-spring  tidal  component  dominate  the  lower  frequency 
band.  The  final  order  of  input  significance  was  surface  elevation,  freshwater 
inflows,  and  wind  components.  A  strong  relationship  was  found  between  sur¬ 
face  elevation  and  wind  components.  Similar  results  also  can  be  obtained  by 
using  cross  spectrum  calculations  with  Parzen  weights.  The  input  significance 
order  is  very  important  when  using  the  correlated  input  decomposition 
method.  Although  the  developed  computer  code  has  that  option,  the  final 
verification  process  used  the  full  matrix  solver. 

In  order  to  conduct  the  model  verification,  several  system  coefficients  were 
determined  after  considering  the  nature  of  the  data.  The  bandwidth  was  the 
most  sensitive  parameter.  This  selection  was  based  on  the  density  of  the  over¬ 
all  spectrum  and  the  accuracy  criteria  of  the  model.  The  final  system 
parameters  were: 

Bandwidth  =  0.00125  cycles'^ 

Maximum  lag  number  =100 

Integration  frequency  range  =  -0.5  to  0.5  cycles"* 

Lag  weight  coefficients  =  1.0 

Memory  lengths  of  response  function  =  -100  to  100  hr 

The  MFRF  which  represent  the  overall  variation  of  salinity  due  to  input 
fluctuations  for  each  point  were  generated  after  verifying  the  model.  The 
multiple  impulse  response  function  (MIRF)  was  obtained  when  the  inverse 
Fourier  transform  of  the  multiple  frequency  response  function  was  conducted. 
The  simulated  salinity  was  generated  by  multiplying  the  MIRF  by  input  series 
over  memory  lengths.  Plate  1  is  a  time  series  comparison  of  the  RMA-10  run 
(solid  line)  and  the  S  run  (dashed  line).  A  350-hr  low-pass  filter  was  used  to 
smooth  the  computed  results.  Aagreement  between  these  two  runs  was  very 
good,  especially  the  low-frequency  variation.  Plate  2  plots  the  corresponding 
residual  series.  The  maximum  computation  error  was  about  0.25  ppt.  The 
final  verification  results  are  shown  as  the  following  four  groups: 

Group  1  :  1990  medium-flow  base  top  layer  (Plates  3-18). 

Group  2  :  1990  medium-flow  base  bottom  layer  (Plates  19-34). 

Group  3  :  1990  medium-flow  phase  I  top  layer  (Plates  35-50). 

Group  4  :  1990  medium-flow  phase  I  bottom  layer  (Plates  51-66). 

The  corresponding  system  simulations  were  summarized  as  another  four 
groups  of  figures  (Plates  67-82,  Plates  83-98,  Plates  99-114,  Plates  115-130). 
The  maximum  variation  of  salinity  regime  during  the  system  simulations  due 
to  the  flow  condition  was  about  2.0  ppt.  No  significant  cyclic  variation  was 
found  in  these  plots.  This  might  be  primarily  attributed  to  using  the  monthly 
mean  input  flow  to  represent  hourly  flow.  This  needs  further  investigation 
with  highly  fluctuating  flow  inputs.  However,  the  advantage  of  system 
simulation  was  time  efficiency.  The  averaged  1-year  hourly  simulation  for 
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each  point  of  interest  only  takes  10  sec  of  Cray  Y-MP  super-computer  central 
processing  unit  time. 
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6  Conclusions 


The  system  modeling  technique  was  used  to  simulate  the  hydrodynamic 
phenomenon  of  salinity  intrusion  in  the  Galveston  Bay  system.  The  modeling 
design  used  numerical  model  results  as  the  output  and  the  same  boundary 
forcing  as  the  numerical  model  for  the  input  to  generate  response  functions. 
While  this  technique  is  limited  to  a  single  point  of  simulation,  it  offers  great 
time-saving  advantages  and  an  easier  model  verification  process.  The 
Galveston  application  showed  very  good  model  verification,  but  the  cyclic 
variation  due  to  flow  condition  change  was  not  significant.  The  success  of 
conducting  the  simulation  using  the  system  model  seems  to  be  highly  cor¬ 
related  to  the  degree  of  dynamic  changes  for  both  the  input  and  output  series. 
Further  research  is  needed  to  identify  the  minimum  sampling  interval  or 
input/output  series  for  the  general  approach. 
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